#Step1: Read the raster
b <- terra::rast("./datal17/temp/cru_ts4.07.1901.2022.tmp.dat.nc")

class(b)

#Step2: Extract time information
time_values <- time(b)

library(terra)
library(sf)
library(stars)
#Step1: Read the raster
b <- rast("./datal17/temp/cru_ts4.07.1901.2022.tmp.dat.nc")

#Step2: Extract time information
time_values <- time(b)

#Step3: Turning strings into time
date_objects <- as.Date(time_values)

#Step4: Select only the dates from the year 1901 and 2022
dates_1901 <- unique(date_objects[format(date_objects, "%Y") == "1901"])
dates_2022 <- unique(date_objects[format(date_objects, "%Y") == "2022"])

#Step5: Find the index of the selected dates in the time_values
selected_indices_1901 <- which(time_values %in% as.Date(dates_1901))
selected_indices_2022 <- which(time_values %in% as.Date(dates_2022))

#Step6: Extract raster values for the selected dates
selected_rasters_1901 <- b[[selected_indices_1901]]
selected_rasters_2022 <- b[[selected_indices_2022]]

#Step7: Extract layer names
layer_names_1901 <- names(selected_rasters_1901)
layer_names_2022 <- names(selected_rasters_2022)

#Step10: Turning the spatraster into a list of raster
stars_list_1901 <- lapply(tmp_layers_1901, function(raster_obj) st_as_stars(raster::raster(raster_obj)))
stars_list_2022 <- lapply(tmp_layers_2022, function(raster_obj) st_as_stars(raster::raster(raster_obj)))

class(tmp_layers_1901)

class(stars_list_1901)

ggplot()+
  geom_stars(data=stars_list_1901[[1]])

ggplot()+
  geom_stars(data=stars_list_1901[[2]])

#Step11: Use lapply to apply st_apply to each stars object to calculate the mean
mean_list_1901 <- lapply(stars_list_1901, function(stars_obj) st_apply(stars_obj, c("x", "y"), mean))
mean_list_2022 <- lapply(stars_list_2022, function(stars_obj) st_apply(stars_obj, c("x", "y"), mean))

#Step12: Combine the list of results into a single stars object
mean_stars_1901 <- do.call(stars::st_as_stars, mean_list_1901)
mean_stars_2022 <- do.call(stars::st_as_stars, mean_list_2022)

#Step8: Find indices of layers that start with "tmp"
#tmp stands for temperature. There are other rasters within these files.
indices_tmp_1901 <- grep("^tmp", layer_names_1901)
tmp_layers_1901 <- selected_rasters_1901[[indices_tmp_1901]]
indices_tmp_2022 <- grep("^tmp", layer_names_2022)
tmp_layers_2022 <- selected_rasters_2022[[indices_tmp_2022]]

#Step9: Changing raster names to temperature
names(tmp_layers_1901)<-dates_1901
names(tmp_layers_2022)<-dates_2022

#Step10: Turning the spatraster into raster
stars_list_1901 <- lapply(tmp_layers_1901, function(obj) st_as_stars(raster::raster(obj)))
stars_list_2022 <- lapply(tmp_layers_2022, function(obj) st_as_stars(raster::raster(obj)))

#Step11: Use lapply to apply st_apply to each stars object to calculate the mean
mean_list_1901 <- lapply(stars_list_1901, function(stars_obj) st_apply(stars_obj, c("x", "y"), mean))
mean_list_2022 <- lapply(stars_list_2022, function(stars_obj) st_apply(stars_obj, c("x", "y"), mean))

#Step12: Combine the list of results into a single stars object
mean_stars_1901 <- do.call(stars::st_as_stars, mean_list_1901)
mean_stars_2022 <- do.call(stars::st_as_stars, mean_list_2022)

#Step13: Changing the name of the rasters
names(mean_stars_1901)<-"avg_tmp_1901"
names(mean_stars_2022)<-"avg_tmp_2022"

library(viridis)
fig<-ggplot()+
  geom_stars(data=mean_stars_1901)+
  scale_fill_viridis(option="turbo", na.value=NA)+
  geom_sf(data=borders, linewidth = 0.01, fill = NA, color = "white", alpha=0.5)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

fig<-ggplot()+
  geom_stars(data=mean_stars_2022)+
  scale_fill_viridis(option="turbo", na.value=NA)+
  geom_sf(data=borders, linewidth = 0.01, fill = NA, color = "white", alpha=0.5)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

library(terra)

# Step 1: Read the raster
b <- rast("./datal17/temp/cru_ts4.07.1901.2022.tmp.dat.nc")

# Step 2: Extract time information
time_values <- time(b)

# Step 3: Turning strings into time
date_objects <- as.Date(time_values)

# Step 4: Get unique years
unique_years <- unique(format(date_objects, "%Y"))
#All years take too long, so we will only choose two years
unique_years <-c("2021", "2022")

# Step 5: Create an empty list to store mean stars for each year
mean_stars_list <- list()

# Step 6: Iterate over each unique year
for (year in unique_years) {
  cat("Calculating average for year", year, "\n")
  # Select dates for the current year
  selected_dates <- unique(date_objects[format(date_objects, "%Y") == year])
  # Find the index of the selected dates in the time_values
  selected_indices <- which(time_values %in% as.Date(selected_dates))
  # Extract raster values for the selected dates
  selected_rasters <- b[[selected_indices]]
  # Extract layer names
  layer_names <- names(selected_rasters)
  # Find indices of layers that start with "tmp"
  indices_tmp <- grep("^tmp", layer_names)
  # Extract temperature layers
  tmp_layers <- selected_rasters[[indices_tmp]]
  # Change raster names to temperature
  names(tmp_layers) <- selected_dates
  # Convert spatraster into raster
  stars_list <- lapply(tmp_layers, function(raster_obj) st_as_stars(raster::raster(raster_obj)))
  # Use lapply to apply st_apply to each stars object to calculate the mean
  mean_list <- lapply(stars_list, function(stars_obj) st_apply(stars_obj, c("x", "y"), mean))
  # Combine the list of results into a single stars object
  mean_stars <- do.call(stars::st_as_stars, mean_list)
  # Change the name of the raster
  names(mean_stars) <- paste("avg_tmp", year, sep = "_")
  # Add the mean stars to the list
  mean_stars_list[[year]] <- mean_stars
  cat("Average calculation for year", year, "completed\n\n")
}

# Combine all the mean stars into a single stars object
# Initialize an empty data frame
result_df <- data.frame()

# Loop through each year and calculate the mean
for (year in names(mean_stars_list)) {
  #Creating an object where we look through the years
  avg_tmp_col_name <- paste0("avg_tmp_", year)
  #Calculating mean and keeping track of the year
  mean_value <- mean(mean_stars_list[[year]][[avg_tmp_col_name]], na.rm = TRUE)
  #Saving everything into a dataframe
  result_df <- rbind(result_df, data.frame(year = as.numeric(year), average_temperature = mean_value))
}

# Writing the df to a dataframe
write.csv(result_df, "./datal17/output/average_temp_2021_2022.csv", row.names = FALSE)

#write.csv(result_df, "./datal17/output/average_temp_1901_2022.csv", row.names = FALSE)
#Reading the file
average_temp<-read.csv("./datal17/output/average_temp_1901_2022.csv")
#Subsetting to below 2001
average_temp_1901_2000<-subset(average_temp, year<2001)
#Calculating the 1901-2000 average
average_temp_1901_2000<-mean(average_temp_1901_2000$average_temperature)
#Calculation the deviations from the 1901-2000 average
average_temp$dev_1901_2000_avg<-average_temp$average_temperature-average_temp_1901_2000

#Mapping the Temperature Deviaitonsover time
ggplot(average_temp, aes(x = year, y = dev_1901_2000_avg)) +
  geom_line() +
  geom_point() +
  theme_bw()+
  labs(title = "Time Trend of Average Temperature",
       x = "Year",
       y = "Average Temperature")

# Create a new column for bar color based on the condition
average_temp$bar_color <- ifelse(average_temp$dev_1901_2000_avg > 0, "red", "blue")

# Create the bar plot
ggplot(average_temp, aes(x = year, y = dev_1901_2000_avg, fill = bar_color)) +
  geom_bar(stat = "identity", color = "black") +
  scale_fill_manual(values = c("blue", "red"), guide="none") +
  theme_bw() +
  labs(title = "Time Trend of Average Temperature",
       x = "Year",
       y = "Difference from 1901-2000 avearge (degrees C)")